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Measuring free energy in parallel tempering Monte Carlo 
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An efficient approach of measuring the absolute free energy in parallel tempering Monte Carlo 
using the exponential averaging method is discussed and the results are compared with those of 
population annealing Monte Carlo using the three-dimensional Edwards-Anderson Ising spin glass 
model as benchmark tests. Numerical results show that parallel tempering, even though uses a 
much less number of temperatures than population annealing, can nevertheless equally efficiently 
measure the absolute free energy by simulating each temperature for longer times. 


I. INTRODUCTION 


Measuring free energy is difficult in Monte Carlo simu¬ 
lations, but nevertheless can be very helpful if available. 
More phenomena and finite size effects can be analyzed 
using free energy. This is not only for the sake of free 
energy, but also as a consequence entropy can be ob¬ 
tained indirectly since energy can be readily measured in 
Monte Carlo simulations. The Bennett acceptance ratio 
method [T] is mainly used to compute the free energy dif¬ 
ference of two different systems but with the same state 
space at the same temperature. To access free energy 
differences of the same system at two different tempera¬ 
tures, exponential averaging (3] and thermodynamic in¬ 
tegration are often used. Thermodynamic integra¬ 
tion has the drawback of integration errors and there¬ 
fore needs data at many temperatures to be accurate. 
When data is available at only few temperatures inter¬ 
polation techniques are often needed. This is especially a 
problem when the specific heat shows complicated multi¬ 
peak structures like in spin glass [5] and protein folding 
[Bj problems due to various reasons like phase transi¬ 
tions and temperature chaos. The exponential averaging 
method however doesn’t have this problem and therefore 
is still beneficial and suitable for such models. 

For systems with complicated free energy landscapes 
like spin glasses, a simple Monte Carlo method is not 
sufficient. Two algorithms that both have a hierarchi¬ 
cal structure are shown to be efficient: parallel temper¬ 
ing pJS] and population annealing [13 El- It has been 
known for more than a decade that population anneal¬ 
ing can accurately estimate the absolute free energy m 
using the exponential averaging method. But slightly 
negative view has been held for parallel tempering due 
to mainly two reasons: the first is that there are too few 
temperatures, therefore the measurement would not be 
very accurate while the other is that the highest tem¬ 
perature is usually finite for parallel tempering while in 
population annealing this can be naturally chosen as in¬ 
finity. The second problem can be trivially solved by 
adding a high temperature stage before the parallel tem¬ 
pering stage using simple single temperature Monte Carlo 
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since thermal equilibration is fast at high temperatures. 
For the first problem, here we show numerically that the 
larger temperature steps in parallel tempering compared 
with population annealing are compensated by simulat¬ 
ing each temperature for longer times, therefore, permit¬ 
ting equally accurate measurements of free energy in par¬ 
allel tempering. Parallel tempering however does suffer 
some minor technical problems, but nevertheless can be 
readily solved to get access to free energy. 

The work is organized as follows: Section |H| introduces 
the Edwards-Anderson model, the free energy estima¬ 
tor using the exponential averaging method and the two- 
stage free energy measurement algorithm in parallel tem¬ 
pering. The results of parallel tempering and population 
annealing are compared in Sec. Ill and the conclusions 
are stated in Sec. Ei 


II. MODELS AND METHODS 

A. The Edwards-Anderson model 

The Edwards-Anderson (EA) Hamiltonian is defined 
as 


T~L ^ ~ Jij SiSj , 
(U) 


(1) 


where Si are the spin degrees of freedom defined on a 
three-dimensional cubic lattice with Sj = ±1. The sum 
over (ij) means sum over the nearest neighbour sites. 
is the coupling between spin s,; and Sj and is indepen¬ 
dently chosen from the standard Gaussian distribution 
with mean zero and variance one. I will refer to each 
disorder realization as a sample. Periodic boundary con¬ 
ditions are used in the simulations of this work. 


B. The free energy estimator 

It has been known that population annealing (PA) can 
naturally estimate the absolute free energy of the EA 
model mm- The method used, exponential averaging, 
is nevertheless general. The ratio of Z(/3 ) and Z(/3), 
which are the partition functions at inverse temperatures 
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where the sum over s is the sum over all states and the 
sum over i is the sum over measured equilibrium states in 
a Monte Carlo run. E is the energy of a micro-state, N 0 
the number of measurements and (...) means a thermal 
average. Take the natural logarithm of both sides 


]nZ(J3')-lnZ(J3)=]nQ(J3,l3') (7) 

-P'F(p') = -PF((3) + lnQ(p,p'), (8) 


where F is the free energy. Note that by sampling equi¬ 
librium states at a temperature, one can integrate the 
free energy between the temperature and a nearby tem¬ 
perature, usually chosen as a lower one. The initial con¬ 
dition of j3F is known at /3 = 0 theoretically. At /? = 0, 
In Z = In Q, where Q is the total number of micro-states 
and Q = 2 n with N = L 3 the total number of spins. If 
we order a set of temperatures as /?o < Pi < ... < fik and 
(3o = 0, then we have 


j=i -1 


- PiF(fii) = > In Q(/3 j+1 ,/3j) + lnfi. (9) 


3=0 


C. Measuring free energy in parallel tempering 

Since parallel tempering (PT) doesn’t usually work 
from /3 = 0, we need to modify the parallel tempering 
algorithm somewhat to get the absolute free energy. One 
simple way is to implement the simulation in two stages 
with a high temperature stage using simple single tem¬ 
perature Monte Carlo and a regular parallel tempering 
Monte Carlo stage. The two stages together will cover 
from fi = 0 to a low temperature of interest. 

• Simple Monte Carlo stage 
Suppose the regular parallel tempering Monte 
Carlo works between T m ; n and T max , the minimum 
and maximum temperatures, respectively. The first 
stage is to use simple single temperature Monte 
Carlo to simulate the system and work from /3 = 0 
to 1/Tmax) not including 1/T max . Then the free 
energy can be integrated from /3 = 0 to 1/T max , 
including 1/T max . 


• Parallel tempering Monte Carlo stage 
Simulate the system using the regular parallel tem¬ 
pering Monte Carlo between T m - m and T max . Then 
the free energy can be further integrated to the low 
temperature spin glass phase. 


In my implementation, I used Np temperatures evenly 
distributed in ft for the first stage while used Nt tem¬ 
peratures evenly distributed in T for the second stage. 
A uniform distribution in /3 for the first stage is easier 
to work with because of the infinite temperature. For 
the second stage, a uniform distribution in T can yield 
better performance of parallel tempering since the swap 
probabilities are more constant at different temperatures. 
Also, no Monte Carlo sweeps were done at j3 = 0, it is 
sufficient to generate random states and make measure¬ 
ments like that of population annealing. In this work, 
the amount of work is counted in terms of Monte Carlo 
sweeps and one Monte Carlo sweep is a sequential up¬ 
date of all the spins for one replica at one temperature 
once. The same number of sweeps were done for both the 
thermal equilibration run and the data collection run. 

It is important to stress here that the only difference of 
this algorithm compared with the regular parallel temper¬ 
ing algorithm is the additional simple Monte Carlo stage. 
Since the free energy landscape is not rough at high tem¬ 
peratures, the autocorrelation time in sweeps therefore 
doesn’t depend on the energy landscape and system size. 
It is not necessary to do as many sweeps as in the second 
stage, where the free energy landscape is rough. There¬ 
fore when the algorithm is well optimized, most of the 
work would be spent in the second stage and the over¬ 
head would be small when adding the first stage. In 
fact, the work spend in the first stage compared with the 
second stage should vanish in the thermodynamic limit. 
In addition, due to the similarity of this algorithm with 
the regular parallel tempering algorithm, it is straight¬ 
forward to modify a regular parallel tempering code to 
collect the free energy data. 

Finally, I want to point out that there is a technical 
issue of overflow when computing Q for large systems if 
Nt is small or there are too many data collection steps 
for parallel tempering. The later problem can be eas¬ 
ily solved since data collection after each Monte Carlo 
sweep is neither necessary nor desired. For the former 
problem, one can subtract a low energy from the energy 
of a state when computing the exponential term and then 
add it back when integrating —f3F. Another method is 
averaging Q on the fly. One can also use reasonably more 
temperatures, but this requires more computational work 
and may also increase the round trip time. Note that 
this problem could appear as well for population anneal¬ 
ing. However, since population annealing naturally uses 
a lot more temperatures and the temperatures are usually 
evenly distributed in /3, no such problem was encountered 
in our studies of the population annealing algorithm at 
least up to size L = 14 down to ff = 3 na. 
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TABLE I: Parameters of the reference runs of PA [14] for dif¬ 
ferent system sizes L with periodic boundary conditions. R 
represents the number of replicas, l//3o is the lowest temper¬ 
ature simulated, Nt the number of temperatures used in the 
annealing schedule, Ns the number of sweeps per temperature 
and M the number of samples. 


L 

R 

1/do 

N t 

N s 

M 

4 

5 10 4 

0.2 

101 

10 

100 

6 

2 10 5 

0.2 

101 

10 

100 

8 

510 s 

0.2 

201 

10 

100 

10 

10 B 

0.2 

301 

10 

100 


TABLE II: Parameters of PT for different system sizes L 
with periodic boundary conditions. Np represents the num¬ 
ber of temperatures in the simple Monte Carlo stage, Nt the 
number of temperatures in the PT stage, T m i n the lowest tem¬ 
perature simulated, T max the highest temperature simulated 
in the PT stage, Ns the number of sweeps per temperature 
and M the number of samples. 


L 

Np 

N t 

T m in 

T max 

N s 

M 

4 

5 

10 

0.2 

2.0 

510 s 

100 

6 

5 

20 

0.2 

2.0 

10® 

100 

8 

10 

30 

0.2 

2.0 

510® 

100 

10 

10 

40 

0.2 

2.0 

5 10 7 

100 


III. RESULTS 

The free energy of each sample was measured using 
both the parallel tempering and population annealing al¬ 
gorithms from infinite temperature to low temperatures 
deep in the spin glass phase. The parallel tempering 
results are compared with the production run of popula¬ 
tion annealing m- The reference simulation parameters 
of population annealing are summarized in Table [I] while 
the simulation parameters of parallel tempering are sum¬ 
marized in Table [TI| Since the free energy results of the 
two algorithms agree to a very high degree of precision, I 
have therefore studied a random subset of samples of Mi 
100 samples per system size. This is sufficient to show 
the equally efficiency of parallel tempering in measuring 
free energy compared with population annealing. In the 
following, I will first make a detailed comparison for a 
single hard sample and then a large scale comparison for 
the two algorithms. 


A. Detailed comparison of a single hard sample 

In this section, I will do a detailed comparison of the 
efficiency of parallel tempering and population annealing 
in measuring the free energy using the hardest sample out 
of about 5000 samples of L = 8 [14]. I will first focus on 
the dimensionless quantity of — (3F at a low temperature 
T = 0.2 and study how the mean and the errorbar of 


the systematic error of —/ 3F change as a function of the 
amount of work W. Then I will study how the relative 
errors of the estimated free energy evolves as a function 
of temperature at a fixed amount of work. The amount of 
work is counted as the total number of sweeps performed 
in the simulation. 

The systematic error A (—/3F) is defined as the dif¬ 
ferences of the estimated (3F and the exact — (3F i.e. 
A (—f3F) = —/ 6F + /3F exact . The exact — f}F is esti¬ 
mated using the reference runs of population annealing 
m- The population annealing comparison data with dif¬ 
ferent amount of work is taken from a previous study of 
finding ground states [T5] . The work is varied by chang¬ 
ing the number of replicas while holding Nt = 101 and 
Ng = 10 constant. The parallel tempering data was done 
using the parameters of Tabel [n] but varying Ns- The 
result is shown in Fig. ]T| One thing we can learn from 
this study is that the accuracy is about the same con¬ 
sidering neither algorithms is very carefully optimized. 
Furthermore, both algorithms underestimate the value 
of — /3F when the amount of work is too small i.e. the 
sample is not in thermal equilibrium. This can be under¬ 
stood as those runs are so small so that some low energy 
states were probably not properly found so that Q was 
underestimated, and therefore also —/3F. 



FIG. 1: Linear-log plot of the systematic error A(—/3F) as a 
function of the amount of work for a hard sample of L = 8 
at T = 0.2. The errorbar is the the standard deviation of 
the —f3F distribution computed from multiple runs, not the 
errorbar of the sample mean of —f3F. The exact value is taken 
from a large scale simulation using population annealing mi¬ 


lt is also important and interesting to study how the 
relative error of the estimated free energy behaves as tem¬ 
perature is lowered due to the integration nature of the 
method. The relative error is defined as minus the ratio of 
the standard deviation ap and mean fip of the estimated 
free energy. The minus sign is to ensure the defined rel¬ 
ative error is positive. The relative error is measured via 
multiple runs in practice. Figure [2] shows this quantity as 
a function of inverse temperature /3 for both PA and PT 
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for the largest run of Fig. [I] Note that there is no trend 
that the relative error grows as temperature is lowered, 
showing the effectiveness of both algorithms. In addi¬ 
tion, the magnitude of the relative errors is also about 
the same. 



P 

FIG. 2: Evolution of the relative error of the estimated free 
energy —uf/^f at a function of the inverse temperature /3 for 
a hard sample of L = 8. The data is the same as the largest 
run of Fig. [I] for both PA and PT. Note that the magnitude 
of the relative error does not grow as temperature is lowered. 


B. A large scale comparison 

Before showing the large scale comparison, I would like 
to show a comparison of the result of — j3F for a typical 
sample of each system size in the whole range of tem¬ 
peratures studied. A typical result is shown in Fig. [3] 
Note that the parallel tempering data falls right on top 
of the population annealing curve, showing the effective¬ 
ness of parallel tempering in measuring free energy in a 
wide range of temperatures. 

To make a comparison of more samples, a scatter plot 
of the free energy per spin / at the lowest simulation 
temperature T = 0.2 of PA and PT is shown in Fig. [I] 
Note that the statistical error compared with the abso¬ 
lute value is too small to be seen in this plot, it is there¬ 
fore interesting to look at the relative error 1— /pt//pa of 
the free energy per spin of parallel tempering against pop¬ 
ulation annealing. The relative error is shown in Fig. [5] 
The errors are well bounded within the accuracy of about 
10 -4 and are well scattered around zero suggesting the 
nature of the errors is essentially statistical and again 
showing PT is as efficient as PA in accurately measuring 
the absolute free energy. 



P 

FIG. 3: Comparison of — j3F for a typical sample of system 
sizes L = 4,6,8 and 10 in a wide range of temperatures. 
The PT data falls right on top of the PA curve, showing the 
effectiveness of PT in measuring free energy. 



-2 - 1.8 - 1.6 - 1 . 4-2 - 1.8 - 1.6 - 1.4 

/pa 

FIG. 4: Scatter plot of the free energy per spin / of PA and 
PT of system sizes L = 4,6, 8 and 10 at T = 0.2. Each point 
represents a sample. 


IV. CONCLUSION 

In this paper, I showed an efficient approach to mea¬ 
sure the absolute free energy in parallel tempering us¬ 
ing the exponential averaging method. The algorithm 
is tested using the three-dimensional EA model and the 
results of parallel tempering agree very well with those 
of population annealing, showing parallel tempering can 
equally efficiently measure free energy as population an¬ 
nealing. There are technical issues with parallel temper¬ 
ing but nevertheless can be readily solved. The fact that 
the reweighting technique works so well suggests that 
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FIG. 5: The same data as in Fig. [4j but plotting the relative 
error 1 — /pt//pa of the free energy per spin of PT against 
PA at T = 0.2. Each point represents a sample. 


many quantities including energy and free energy can be 
accurately estimated using the reweighting technique at 
many temperatures with little overhead in parallel tem¬ 
pering without using interpolation techniques. Finally, 
I hope that this simple and direct access to free energy 
in parallel tempering can make the spin glass and other 
relevant research fields potentially richer. 
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